SUBROUTINE olsreg(vsize,Nreg,yin,xin,betaout,Vresout,E2out,out_err)
  USE Commonvars

  IMPLICIT NONE
  INCLUDE 'mpif.h'
  
  INTEGER, INTENT(IN)  :: vsize, Nreg  
  REAL(8), INTENT(IN)  :: yin(vsize), xin(vsize,Nreg)
  REAL(8), INTENT(OUT) :: Vresout,E2out,betaout(Nreg),out_err
  REAL(8), ALLOCATABLE :: res(:), work(:)
  REAL(8)              :: Mres, xpx(Nreg,Nreg), xpy(Nreg), invxx(Nreg,Nreg)
  INTEGER              :: info, ii, jj, lwork
  CHARACTER            :: uplo='U'

  out_err = 0d0

  ALLOCATE(res(vsize))
  res = 0.0d0
  
  xpx = 0.0d0
  xpy = 0.0d0
  
  xpx = MATMUL(TRANSPOSE(xin(:,:)),xin(:,:))/DBLE(SIZE(yin))
  xpy = MATMUL(TRANSPOSE(xin(:,:)),yin(:))/DBLE(SIZE(yin))
  
  invxx=0.0d0
  DO jj = 1, Nreg
     invxx(jj,jj)=1.0d0
  END DO
  
  LWORK = Nreg
  IF (ALLOCATED(WORK)) DEALLOCATE(WORK)
  ALLOCATE(WORK(LWORK))
  CALL DPOSV(uplo,Nreg,Nreg,xpx,Nreg,invxx,Nreg,info)
  
  IF (info /= 0) THEN
     DO jj = 1, Nreg
        IF (myrank == 0) WRITE(*,*) 'xpx', xpx(jj,:)
     ENDDO
     WRITE(*,*) 'MATRIX COULD NOT INVERT simulatemun', 'info', info, 'myrank', myrank, 'Nreg', Nreg
     DO ii = 1, vsize
        IF (myrank == 0) WRITE(*,*) 'y', yin(ii), 'x', xin(ii,:)
     END DO

     out_err = 1d0

  END IF
  
  betaout = MATMUL(invxx,xpy)

  res = yin(:) - MATMUL(xin(:,:),betaout(:))
  
!!$  Rsquared = (SUM((yin(:)-SUM(yin(:))/SIZE(yin(:)))**2d0) - SUM(res(:)**2d0)) / SUM((yin(:)-SUM(yin(:))/SIZE(yin(:)))**2d0)
!!$
!!$  if (flag_here==6) write(*,*) 'Rsquared', Rsquared

  Mres = SUM(res)/DBLE(vsize)
  Vresout = SQRT(SUM((res(:)-Mres)**2d0)/DBLE(vsize))
  E2out = SUM(res(:)**2d0)/DBLE(vsize)
  
  DEALLOCATE(res)
  
END SUBROUTINE olsreg
